#! /usr/bin/env python
# -*- coding: utf-8 -*-

"""
Расчет грида плотности распределения точек
Под точками подразумевается нормали к плоскостям трещин
Вариант для прямоуголной диаграммы
"""
from math import sqrt, exp, pi, cos, radians
from plane import Plane
from grid_plane import PlaneGrid
import sys
from optparse import OptionParser

def gauss(m):
    # return value of a gauss distribution.
    # ( 1 / (sigma*sqrt(2*pi)) ) * exp (-1* (m)**2/(2*sigma**2) )
    # it use following constants:
    # sigma_pi = 1 / (sigma*sqrt(2*pi))
    # sigma_2 = 2*sigma**2
    #return sigma_pi * exp ( -1*(m)**2/sigma_2)
    return exp ( -1*(m)**2/sigma_2)
    #if m>sgm1: return 0.0
    #return (sgm1-m)/sgm1
    #if m>sgm1: return 0.0
    #return 1.0

def segment_square(dip,incr):
    return 2*pi*1.0*(cos(radians(dip)) - cos(radians(dip+incr))) * (incr*1.0/360.0)


def line_parse(line):
    line = line.expandtabs(1)
    line = [x for x in line.split(' ') if x != '' ]    
    DIP, DIR = 0, 0    
    
    try:
        DIR, DIP = int(line.pop(0)), int(line.pop(0))
    except:
        #sys.stdout.write('WARNING! Cannot parse string: %s\n' % line)
        return None   
    line = Plane(DIR, DIP)
    return line

# Parsing of an options
optparser= OptionParser()
optparser.add_option('-s', '--sigma', dest='sigma',
                     default=10)
optparser.add_option('-i', '--increment', dest='increment',
                     default=5)
optparser.add_option('-j', '--intensity', dest='intensity',
                     action='store_true',
                     default=False)
options, args = optparser.parse_args(sys.argv[1:])

lines = sys.stdin.read().split('\n')

# parsing of the data
data=[]
if options.intensity:
    # incoming data are planes with intensity
    for line in lines:
        vals = [float(x) for x in line.split(" ") if x !='']
        if len(vals) !=0:
            data.append([Plane(vals[0],vals[1]), vals[2]])
else:
    for line in lines:
        val = line_parse(line)
        if val != None: data.append(val)  

sigma = float(options.sigma)
increment = int(options.increment)
dirsize = 360.0/increment # number of raster
dipsize = 90.0/increment

grid = PlaneGrid(incr=[increment, increment])

# constants for calculating gauss
global sigma_pi
global sigma_2
global sgm
global sgm1
sgm1 = sigma*1.0
sgm = sigma*2
#sigma_pi = 1 / (sigma*sqrt(2*pi))
sigma_2 = 2*sigma**2
i = .0
max_i = len(data)*(360.0/increment)*(90.0/increment)

# calculate fuction value at a zero-point to normalize all the other values
#gauss0 = 1/gauss(0,sigma)


mval=0.0
if options.intensity:
    for ddir in range(0, 360, increment):
        for ddip in range(0, 90, increment):
            for vals in data:
                plane=vals[0]
                intensity=vals[1]
                i+=1
                sys.stderr.write('%g\r' % (i/max_i))
                cell_plane=Plane(ddir, ddip)            
                angle = cell_plane.normal().return_angle_between(plane.normal())
                if angle>90: angle = 180 - angle
                grid.add(plane=cell_plane, val=gauss(angle)*intensity)
                mval = max(grid.get(plane=cell_plane),mval)
else:
    for ddir in range(0, 360, increment):
        for ddip in range(0, 90, increment):
            for plane in data:
                i+=1
                sys.stderr.write('%g\r' % (i/max_i))
                cell_plane=Plane(ddir, ddip)            
                angle = cell_plane.normal().return_angle_between(plane.normal())
                if angle>90: angle = 180 - angle
                grid.add(plane=cell_plane, val=gauss(angle))
                mval = max(grid.get(plane=cell_plane),mval)

sys.stdout.write('%2.2f\n' % (increment))
sys.stdout.write('%2.2f\n' % (mval))
for ddir in range(0, 360, increment):
    for ddip in range(0, 90, increment):
        a = Plane(ddir,ddip)
        sys.stdout.write('%d\t%d\t%g\n' % (a.dir, a.dip, grid.get(plane=a)))
